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ABSTRACT 

The inspiral and merger of binary black holes will likely involve black holes with both unequal masses 
and arbitrary spins. The gravitational radiation emitted by these binaries will carry angular as well 
as linear momentum. A net flux of emitted linear momentum implies that the black hole produced by 
the merger will experience a recoil or kick. Previous studies have focused on the recoil velocity from 
unequal mass, non-spinning binaries. We present results from simulations of equal mass but spinning 
black hole binaries and show how a significant gravitational recoil can also be obtained in these 
situations. We consider the case of black holes with opposite spins of magnitude a aligned/anti-aligned 
with the orbital angular momentum, with a the dimensionless spin parameters of the individual holes. 
For the initial setups under consideration, we find a recoil velocity of V — 475 km s -1 a . Supermassive 
black hole mergers producing kicks of this magnitude could result in the ejection from the cores of 
dwarf galaxies of the final hole produced by the collision. 

Subject headings: black hole physics — galaxies: nuclei — gravitation — gravitational waves - 
relativity 
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1. INTRODUCTION 

There is ample observational evidence that super- 
massive black h oles (SMBHs) are common at the cen- 
ters o f galaxies (|Richstone et al.l 119981 : iMagorrian et al.l 
119981) . with masses in the range 10 5 - 1O 9 M . These 
SMBHs are involved in exciting astrophysical phenom- 
ena. For instance, there is a remarkable, not com- 
pletely understood, correlation between the velocity dis- 
persion of the bulge of the host galaxy and the mass 
of the SMBH (jFerrarese fc Merritu 120001 ). There is also 
indication of a correlation of the mass of t he SMBH 
with the mass of the host dark matter halo (|Ferraresd 
2002). An interesting aspect of SMBH growth arises 
as a consequence of hierarchical cold dark matter cos- 
mologies, in which large-scale structures are formed by 
mergers. SMBHs would then grow both by gas accretion 
and by coalescence with other SM BHs (brought together 
when their host gala xies collide (|Volonteri et al.ll2003t 
iBegelman et al.lll980[ )L The work in this paper focuses 
on one aspect of the merger of SMBHs, the kick in the 
final SMBH. 

The late inspiral and merger of SMBHs produces 
extremely energetic gravitational radiation, which will 
be observable by the planned spa c e-bas e d gravitationa l 
wave antenna LISA (|Danzmannl 120031: iPrincd I2003D . 
Gravitational radiation produced during the inspiral and 
merger of black holes (BHs) not only carries energy 
with it, but, except in special-symmetry cases, can also 
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transport net linear and angular momentum. For in- 
stance, in the merger of unequal mass SMBHs, a net 
flux of linear m o mentum will be emitted by the sys- 
tem (|Pereslll962t lBekensteinlll973l ). As a consequence, 
the final BH will experience a gravitational recoil or kick. 
There are observations that hint at such scenarios, in 
which a SMBH has been eje cted in an ongoing galaxy 
merger (jHaehnelt et al.l[2006T l . (An alternative explana- 
tion could be that the ejection is due to gravitational 
slingshot of three or more SMBHs in the merger.) It is 
then very important to get good estimates of recoil veloc- 
ities in BH mergers. These estimates have a profound ef- 
fect on the understanding of the demogra phics of SMBHs 
at the cores of galaxies , their growth ( Haima n1 [200l 
and their merger rates (|M 1C1C et all 12006|). Knowledge 
of the conditions under which kicks are produced could 
also help explain the absence of massive BHs in dwarf 
galaxies and stella r clusters (jMadau fc Quataertl 120041 : 
iMerritt et al.l l2004) , and could determine the population 
of BHs in the interstellar and intergalactic medium. 

Gravitational recoil estimates of unequal mass bina- 
ries have been addressed using both analytic and full nu- 
merical relativity a pproaches. The first quasi-Newtonian 
analy tic studies (|FitchettJ |1983t iFitchett k, Detweilerl 
1984) pro d uced k ick velocities as large as < 



1500 km s~ 



Wisemanl (Il992f) and more recent l y iBlanchet et al.l 
(I2005h and iDamouT &: Gopakumarl (|2006t ) improved 
these estimates by including post-Newtonian (PN) ef- 
fects. The maximum kick in these studies was found to 
be in the range of ~ 74 — 250 km s , and it occurred 
for i] ~ 0.2, where r\ = M 1 M 2 /(M 1 + M 2 ) 2 is the sym- 
metrized mass ratio parameter. (This corresponds to a 
mass ratio q = M\fMi ~ 0.38.) These analytic PN stud- 
ies also showed that the final value of the kick is mostly 
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accumulated during the merger or plunge phase of the 
binary. Since the plunge phase is beyond the limit of ap- 
plicability of PN approximations , the results can only be 
taken as "best-bet estimates" (iDamour fc Gopakumarl 
[2001 . 

There are two semi-analytic studies that in principle 
had a better handle on the plunge phase. ICampanellil 
(2005) obtained kick velocities of ~ 300 km s" 1 using 
the Lazarus approach, a framework (B aker et al.l 120021) 
that combines full numerical relativity an d close-limit ap- 
proxi mation (CLA) pertur bation theory (IPrice fc Pullinl 
19941). More recen tly, ISopuerta et all (|2007| ) and 



Sopuerta et alj (|2006h combined PN estimates during 
the inspiral with kick estimates using the CLA. The 
maximum recoil obtained in this work was ~ 167(1 + 
e)km s~ , with e the eccentricity of the binary. Fi- 
nally, full n umerical relativity studi es have also be e n car - 
ried out by iHerrmann et alt (|2006f) . iBaker et all (|2006l ) 
and iGonzalez et alj (|2006f h Only full numerical rela- 
tivity approaches provide accurate estimates of kicks 
since they correctly handle the non-linear behavior of the 
plunge. The most com prehensive study so far is that by 
IGonzalez et al.1 (|2006f ) , in which a maximum kick velocity 
of ~ 175 km s~ was obtained also for r\ ~ 0.2 (g ~ 0.38), 
consistent with PN studies. What i s interesting is that 
the findings of ISopuerta et a l. (2007) based on the CLA 
are rem arkably close to the fu ll numerical relativity re- 
sults by I Gonzalez et ahl (|2006f ). supporting the view that 
the kick is mostly due to the linear momentum emit- 
ted during the plunge, where the CLA ha s been demon- 
strat ed to provide a good approximation (|Anninos et al.l 
fl995h . To our knowledg e, the only kick stu dy involving 
spinning BHs is that by iFavata et alj (|2004D . They con- 
sidered the case of an extreme-mass-ratio system with a 
spinning SMBH. Using BH perturbation theory they es- 
timated kick velocities of ~ 100 — 200 km s _1 . Head-on 
collision s of spinnin g BHs have also been recently con- 
sidered (|Choill2007h . 

To help us understand our computational results, we 
present next a rough order-of-magnitude estimate of the 
kicks one should expect. Note first that there must be 
some asymmetry between the BHs in order for there to 
be asymmetric radiation which can lead to kicks. Thus, 
in the non-spinning case, unequal masses are required; 
here we consider binaries of equal masses, but different 
spin (magnitude or direction) . The kick is expected to in- 
crease as the relevant spin increases, but especially sym- 
metric cases will still show zero kick (e.g. when the BHs 
have their spins aligned parallel to the orbital angular 
momentum). The order-of-magnitude estimate can be 
obtai ned from the r adiative linear momentum loss for- 
mula (|ThorneHl98d lKidderll!995ah . Excluding non-spm 
terms, this formula reads 
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Here Jjj and 1^ are respectively the mass quadrupole 
and octupole. Similarly, Hij and Hijk are the spin 
quadrupole and octupole, respectively. In |T]), a super- 
index ( n > denotes an nth-time derivative. Clearly equa- 
tion Jl]) predicts a periodic force for exactly circular or- 



bits. As the BHs spiral together the strength of the peri- 
odic kick increases, so we estimate the kick from the last 
half orbit before merger. 

Consider a binary system consisting of BHs in circu- 
lar orbit with equal masses (Mi = A/2 = M/2). In the 
absence of spin this would produce no kick, but here we 
set data with each BH having spin perpendicular to the 
orbit, the spins oppositely directed, each with dimen- 
sionless Kerr spin parameter a (0 < a < 1). This is the 
configuration we use below for our computational evalua- 
tion of the kick. The calculation of the mass quadrupole 
is familiar, and for circular orbits in the xy plane with 
orbital angular velocity lo and coordinate separation d 
gives nonzero values: 



li a J=2M<F lj 3 sin(2wt) 
I^=-2M d 2 u 3 cos(2urt) 
2M d 2 lo 3 sin (2wt) . 



r(3) 
yy 



(2) 



The spin quadrupole can be most easily calculated 
by imagining a spin dipole (charges ±M/2, separation 
a M/2) and conceptually taking the limit at the end. The 
result is 



Uf) = - M 2 d a lo 3 sin (ut) 
H$ = --M 2 dauj 3 cos (cut) . 



(3) 



Inserting and into the first term in equation JT]) 
gives 

dP x 



M s d 6 auj b sin (ut) 



dt 45" 

dpy 8 , ,0 fi 

= M i d i au b cos (ut) . 

dt 45 v ' 



(4) 



Notice that the force is in the plane of the orbit and 
rotates with the orbit. The average over half a cycle 
is 2/tt, so equation (J5J) is a good estimate for any half 
cycle as the orbit spirals in. The total force can then be 
approximated as 

dP 16 



-M 3 d 3 auj 6 . 
dt 45 7T 

Compare this to the total luminosity: 
dE 2 



dt 



-M 2 d 4 LU 6 . 



(5) 



(6) 



Thus the asymmetry in radiation that contributes to the 
kick is 



dP _ dP ,dE _ 4a M 
dE ~ ~dt'~dt~9^~d 



(7) 



which is ~ 0.02 for dimensionless spin parameter a ~ 1/2 
and d ~ 6M. (The latter is an estimate of the separation 
near the "last orbit" .) If AE is the total energy radiated 
by the binary, an estimate of the (half orbit) kick is 

AE\ 



V > 



dP_ 

dE 



For another estimate, we note that IFavata et al 
(2004) specialized the PN equation (3.31) in iKidder 
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(|1995bh to the case of circular orbit with spins paral- 
lel and anti-parallel to the orbital angular momentum. 
The resultant kick velocity is given by 

V = V q + 88^s-i(f s f> al >^)(™). (9) 

\ J SO. max. / \ ^term / 

Above V q is the contribution to the kick that de- 
pends only on the mass ratio q; this contribution van- 
ishes for equal mass binaries (q = 1). The radius 
'"term is the separation at which gravitational radiation 
terminates. The scaling function in equation ^ is 
given by fso(q, Oi, a 2 ) = q 2 {a 2 - £?Oi)/(l + q) 5 with 
/so.max = /so(l,±l,Tl) = 1/16. Therefore, for 
the cases we have investigated, equation © reduces to 
V = 883 km s~ a (2 M/r tcrm ), comparable to our esti- 
mate above, for reasonable choices of rt erm . 

There is another effect similar to the pulsar kick mech- 
anism described bv lHarrison fe T adcmarul (|1975f) . It in- 
volves explicit retardation effects (so is not captured in 
the multipole expression of equation ([1])), and gives es- 
timates of similarly sized kicks. We shall see that full 
numerical relativity simulations give comparable kicks to 
this estimate. 

The paper is organized as follows: In SecEJ we present 
the computational methodology and details of how the 
initial data were constructed. Sec. [3] gives details of the 
method to estimate kicks. Code tests and a convergence 
analysis are given in Sec.[U The gravitational recoil esti- 
mates are presented in Sec. [5l with conclusions given in 
Sec. El 

2. COMPUTATIONAL METHODOLOGY AND INITIAL 
DATA 

The numerical simulations of binary black holes (BBH) 
in our work were obtained following the Moving Punc- 
ture Recipe (MPR). The essence of this recipe is: (A) 
a particular formulation of the Einstein field equations 
and (B) a set of coordinate or gauge conditions for up- 
dating field variables during evolution as well as for han- 
dling the BH singularities. The form of the evolutions 
required by the MPR is the so- called BSSN 3+1 for- 
mulation of Einstein's equations dNakamura et al.lll987t 
IShibata fe Nakamur a 1995; Baumgarte fc Shapirdll999f) . 
A derivation of the BSSN equations and a few exam- 
ple s of their applications can be found in the review 
bv lBaumgarte fc Shapirol (|2003fh 

In addition to the form of the evolution equations, the 
success of t he MPR is due to the coordinate o r gauge 
conditions (lAlcubierre et all 120031 : iBaker et all l2006bl: 
ICampanelli et alj l2006al ) . The MPR gauge conditions 
are equations that determine the lapse function a and the 
shift vector (3 l . The lapse is a "local" measure of proper 
time, and the shift vector encap sulates the freedom of 
label ing events at a given time (jBaumgarte fc Shapirol 
2003). The explicit form of the evolution equations for 
the lapse and shift in the MPR are d^a = —2aK and for 
the shift d (3\ = 3/45* and d B i = fyF - £B\ where 
do = dt — ftdj. K is the trace of the extrinsic cur- 
vature, T l is the trace of the conformal connection and 
£ = 2 is a free dissipative parameter. The importance of 
these gauge conditions is twofold: first, they avoid the 
need of excising the BH singularity from the computa- 
tional domain since they effectively halt the evolution 



TABLE 1 
Initial Data Parameters 



Model 


x/M 


P/M 


S/M 2 


mi/M 


m 2 /M 


E/M 


S0.05 
S0.10 
S0.15 
S0.20 


2.95 
2.98 
3.05 
3.15 


0.13983 
0.13842 
0.13547 
0.13095 


0.05 
0.10 
0.15 
0.20 


0.4683 
0.4436 
0.3951 
0.2968 


0.4685 
0.4438 
0.3953 
0.2970 


0.98445 
0.98455 
0.98473 
0.98499 



(i.e. the lapse function a va nishes) near the BH sin- 
gularity (Han nam et al.ll2006f ). Second, they allow for 
movement of the BH or puncture throughout the compu- 
tational domain while freezing the evolution inside the 
BH horizon. 

The code used for this work was produced b y the 
Kranc code generation package dHusa et al.l [2006f ) . the 
Cactus infra structure (|Cactusl 120071) for parallelization 
and Carpet (|Schnetter et al.ll2004h for mesh refinement. 
The code is based on fourth order accurate finite differ- 
encing of spatial operators and uses 4th order Runge- 
Kutta for time integration with a Courant factor of 0.5. 

Th e initial data use punctures (|Brandt fc Brugmannl 
Il997f ) to represent BHs. In Einstein's theory, initial data 
are not completely freely specifiable; they must satisfy 
the Hamiltonian and momentu m constraints . We use 
the spectral code developed by lAnsorg et al.l (|2004h to 
solve these constraints. The initial free-data (e.g. an- 
gular momentum, spins, masses, separations ) are chosen 
according to the effective potential method (|CookllT994l : 
lBaumgartell2000D . This method yields BBH initial data 
sets representing BBHs in quasi-circular orbit. In general 
terms, the effective potential method consists of minimiz- 
ing the "binding energy" of the binary to determine the 
BBH parameters. 

Table Q] contains the BBH parameters of our simula- 
tions. The BHs are located at positions (ztx/M, 0, 0), 
have linear momentum (+P/M, 0, 0), spin (0, 0, ±S/M 2 ) 
and bare puncture masses toi^/M, with M = M\ + M 2 
the total mass of the binary. Notice that the bare punc- 
ture masses are slightly different. The reason for this 
difference is because of the spin contribution to the mass 
of each hole (measured from the area of their appar- 
ent horizons); in order to keep the individual masses 
of the BHs, M\ and M 2 , equal, (slight) adjustments to 
the bare masses are necessary. The configurations are 
such that the total angular momentum is for all cases 
J/(/xM) = 3.3 with n = M 1 M 2 /(M 1 + M 2 ). It is im- 
portant to notice that S is not the Kerr spin parameter 
< CLKerr < TnBH typically associated with rotating 
BHs. The dimensionless spin parameter for each BH is 
given by a± t 2 = S/M 2 2 with Mi 2 = M/2. The cases con- 
sidered here, S/M 2 = {0.05,0.10,0.15,0.20}, correspond 
to a = {0.2,0.4,0.6,0.8}, respectively. For reference, the 
total ADM mass E/M in the initial data is also reported 
in Table □ 

The computational grids consist of a nested set of 10 
refinement levels, with the finest mesh having resolution 
h = M/40. This resolution translates into a resolution 
of about h = m/19 — m/12, with respect to the bare 
mass m of the puncture according to Table [TJ The min- 
imal resoluti on found to be adequate for spinning cases 
according to ICampanelli et alj (|2006bD is ft. < M/30. In 
our h = M/40 simulations there are 4 refinement levels 
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of 58 3 grid-points nested within 6 levels of 102 3 grid- 
points. During the evolution the shape and number of 
grid-points per refinement level vary as the centers of the 
grids track the positions of the black holes. The coarsest 
mesh is kept fixed and extends to 650M from the origin 
in each direction. 

3. GRAVITATIONAL RECOIL 

The gravitational recoil is computed from the rate of 
change of linear momentum 



dP i 
~~dt 



= lim 



d 2 E 

7 

dVldt 



r 2 dQ 



(10) 



which is determined by the fluxes of energy E and lin- 
ear momentum P L (n l is the unit normal to the sphere). 
In order to compute the recoil velocity, the Newtonian 
momentum relation is used, V % — P l /M. 

In terms of ^4, the component of the Weyl curvature 
tensor representing outgoin g radiation, equation (I10| 
reads ([Newman fc Todlll980t ) 



dP l 
~dT 



lim ■ 

r — >oo 



[4tt/ 


j * 4 dt' 


n l r 2 dftj 




J —OO 





(11) 



Equation (fTTj) is applied at a finite radius r > 30M 
away from the "center of mass" of the binary but far 
enough from the boundary of the computational domain 
to avoid t he effects from spurio us reflection from the 
boundary (|Zlochower et al.l l2Q05h . The Weyl scalar * 4 
is computed in the bulk of the computational domain 
and is then projected onto the sphere and used in the 
computation of equation (jll|) . 

We also estimate the gravitational recoil using a mode 
decomposition. Instead of constructing 'J' 4 in the bulk 
of the computational domain and interpolating it on a 
sphere to be used in equation pip , we decompose ^4 into 
spin- weight —2 spherical harmonics and then compute 
the recoil. That is, one first constructs the coefficients 
-2Ctm such that 



*4 = X! -2Ct m {t,r) - 2 Yt m (6,tp). 



(12) 



Given these coefficients, the gravitational recoil is given 

by 

dP i 

' =^2(£,m\lm) (13) 

£m£ni 



dt 



where (£, m\£, m) represents the contribution to dP l /dt 
from the overlap 



(£, m\£, m) oc Re 



-iYl m -iY lih d£l 



(14) 

with _ 2 C£m = /-oo -^Pi m dt' . This mode-overlap de- 
composition has the advantage that the contribution 
from different overlapping modes can be studied indi- 
vidually. 

There is an important issue to keep in mind when us- 
ing both equations (fTTj) and (fT3j) to estimate kicks. It 
is well known that initial data in BBH simulations con- 
tain spurious radiation. Fortunately, this radiation does 
not seem to have a significant effect on the dynamics of 
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Fig. 1. — Fluxes of energy dE/dt, linear momentum dP l /dt 
and angular momentum dj/dt as a function of time for the SO. 10 
(a = 0.4) case. The vertical line at 60 M denotes t m i n , the lower 
limit of the time integration used to estimates kicks which avoids 
contamination from the spurious radiation in the initial data. 

the binary. However, because of the time-integration in- 
volved in the kick formulas, the estimates are affected 
by the spurious radiation. To alleviate this problem, we 
set the lower limit in the time integral to be i m ; n and 
choose /j m j n as the time after which the spurious burst 
has passed. As an example, Figure [TJ displays the fluxes 
of energy dE/dt, linear momentum dP % j dt and angular 
momentum dJ/dt through the detector at rdet = 30 M 
for the ST). 10 case. It is clear from these rates that there 
is a spurious burst from the initial data for t < 50 M. In 
particular, notice the effect on dP % jdt at early times. The 
line at i m i n = 60 M shows our choice for this cut-off. The 
precise choice of 7j m j n is not important, as long as the ini- 
tial spurious burst is eliminated and i m i n is not too close 
to the time when the amplitude of the gravitational wave 
becomes relevant. Since we use several locations ("detec- 
tors") at different radii to compute fluxes, the value of 
tmin is adjusted as t mm = 30 M + rdet, where r^et denotes 
the detector radius. Note the smallness of dP z /dt from 
FigureQ] It translates to velocities of ~ 0.2 km s _1 ; thus, 
we will not plot V z in subsequent figures. 

Another important check when computing kicks us- 
ing equations (TTT|) and ([13"]) is the dependence of the re- 
sults on the extraction radius r^ ct . The kick formulas 
are in principle valid in the limit r — > 00, but one applies 
them at a finite extraction radius rdet where there is suf- 
ficient resolution. Figure [2] shows the recoil velocity as 
a function of time computed at different detector radii, 
Tiet/M — (30, 40, 50). The time dependence of the veloc- 
ities has not been adjusted by the lag in arrival times at 
each detector. Although small, one can see from Figure^] 
that there a is slight sensitivity of the extracted kick ve- 
locity to the location of the detector for the ranges we 
considered. This variation is within the error estimates 
of our kicks. The origin of this dependence of the ex- 
tracted kick on the detector location could be numerical 
(e.g. outer boundary, mesh refinement interfaces, etc.) or 
due to the redshift and tail effects. 3 

4. CODE TESTS AND WAVEFORM CONVERGENCE 

We have tested that our code produces a sufficient level 
of convergence for equal mass, non-spinning BH binaries 
that we are confident in the res ults. In particular, we 
have carried out extensive tests (Shoema ker et al1l2007f ) 

3 We thank the anonymous referee for bringing this to our at- 
tention. 
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150 
time (M) 

Fig. 2. — Recoil velocity V x and V y computed from different de- 
tector locations for 50.10 with resolution h = M/40. The detectors 
were located at r dBt /M = (30, 40, 50). 



for the Rl run in Bak er et alJ (|2006af ) and found res- 
olution ranges that yield between 3rd- and 4th- order 
convergence. Also as a code test, we carried out a non- 
spinning, unequal mass simulation for r; = 0.23. The kick 
obtained from t his ru n (~ 130km s _1 ) matches that by 
iGonzalez et alJ (|2006l ) . Because the BBH setups in our 
present work have no symmetries, the computational cost 
of each simulation is high (for our h = M/40 resolution 
runs the cost is ~ 44 hours on 32 CPU cores for a total 
of about ~ 1400 CPU hours on a supercomputer), so to 
demonstrate convergence our runs were limited to reso- 
lutions h < M/40. We present convergence results for 
the 50.10 case; the other cases have similar behavior. 

Figure [3] shows the amplitude of the dominant £ = 
2, to = 2 mode of ^4. The top panel of the fig- 
ure displays the mode at the three different resolutions 
(h/M = 1/32, 1/35, 1/40). The bottom panel shows the 
coarse-medium ("c-m") differences and the medium- fine 
("m-f") differences reseated for 2nd, 3rd and 4th order. 
As the plot shows, this mode converges between 3rd and 
4th order. In our convergence studies for other systems 
(e.g. equal mass BHs) getting closer to 4th-order con- 
vergence required at least a factor of two between the 
coarsest and finest resolution. Given the range of reso- 
lutions that we are able to do for the present study, the 
deterioration of our convergence should not be surprising. 
Nonetheless, we believe that the observed level of conver- 
gence in our simulations will not affect the astrophysical 
implications of the magnitude of our kick estimates. 

As a check of our implementation of the kick extrac- 
tion, Figure |4] compares the recoil velocity computed 
from equation (|lip and equation (fT3"|) for the case 50.10 
with resolution h = M/40. For equation (| 13[) . we include 
up to I = 4 modes. It is evident from this plot that with 
the modes I < 4 one can reconstruct most of the total 
recoil velocity. 

5. RESULTS 

First, we present the main results of our work, namely 
the kick estimates together with the radiated energy and 
angular momentum, followed by a discussion of conver- 
gence and a mode analysis of the kicks. 

5.1. Kicks and Radiated Energy and Momentum 

The core results of our work are summarized in Ta- 
ble El Table [2] lists the values for the total recoil V, en- 
ergy AE and angular momentum A J radiated for each of 
the cases considered. The reported values were obtained 
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Fig. 3. — The amplitude of the dominant I = 2, m = 2 mode 
of \p4 for the case SO. 10 (a = 0.4). The top plot shows the mode 
at three different resolutions (h/M = 1/32,1/35,1/40), while the 
bottom shows the small differences between the medium-coarse 
("c-m") and the medium-fine ("m-f") simulations rescaled for 2nd, 
3rd and 4th order. The waveform is between 3rd- and 4th-order 
convergent. 
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Fig. 4. — R ecoil velocity V x a nd V y versus time computed from 
equation JTTJ and equation (1 1 3 I t for the SO. 10 model with resolu- 
tion h = M/A0 extracted at r det = 30 M. V z is below 0.2 km/s 
and hence is not shown. The insets labeled "dif ferences" show the 
difference between the recoil from equation mil l and equation (1 1 3 P 
with modes up to and including I = 4 



TABLE 2 
Radiated Quantities 



Model 


a 


V(km s" 1 ) 


AE(%) 


AJ(%) 


S0.05 


0.2 


96 ± 7 


3.24 


26.82 


S0.10 


0.1 


190 ± 10 


3.30 


27.05 


S0.15 


0.6 


285 ± 12 


3.33 


27.12 


S0.20 


0.8 


392 ± 33 


3.34 


26.83 



with resolutions h — M/40 and extracted at rdet = 40 M. 
For reference, we include also the dimensionless spin pa- 
rameter a. Figure [5] displays the recoil velocity V as a 
function of the dimensionless spin parameter a for all the 
resolutions used in our simulations. Solid circles denote 
resolutions h = M/40, diamonds resolutions h = M/32 
and inverted triangles resolutions h = M/30. The error 
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Fig. 5. — Magnitude of the recoil velocity V as a function of the 
dimensionless spin parameter a. Solid circles are for resolutions h = 
M/40, diamonds for resolutions h = M/32 and inverted triangles 
for resolutions h = M/30. In each case, the results at different 
resolution cluster more tightly than the conservatively estimated 
error bars (Table lo- 
bars correspond to the conservatively estimated errors 
listed in Table [21 and are larger than the actual scatter 
of the results at different resolution. 

In order to estimate these errors, for each spin case, 
we perform Richardson error estimates of the total recoil 
velocity V assuming 2nd order convergence. We then in- 
crease these errors to take into account factors such as the 
deterioration of convergence in the weak mode-overlaps 
(see below). We believe these are conservative best-guess 
errors that could be reduced with, among other things, 
higher resolution. 

Note in Figure [5] the linear dependence of the magni- 
tude of the kick velocity V on the spin parameter, as 
expected from the multipole example in Section [TJ A fit 
to the data yields V — 475 km s _1 a . 

An interesting aspect of the spin configuration we have 
considered is the fraction of radiated energy AE and an- 
gular momentum A J. The fraction radiated is approxi- 
mately constant within the accuracy of our simulations. 
One possible reason why AE and AJ do not seem to 
depend on the spins of the holes could be due to the set 
up of our initial data. By construction, the four cases 
we considered have the same total initial angular mo- 
mentum J//zM = 3.3. In our case with spins oppositely 
directed and with equal magnitude the variations in the 
total ADM energy are < 0.05%, as can be seen from 
Table [U 

5.2. Mode Analysis and Convergence 
With the kick formula (fT3"]l . we were able to investi- 
gate the contribution of each mode-overlap (£, m\£, fh) to 
the total recoil velocity. Figure [6] shows the contribution 
that each mode-overlap makes to the total kick veloc- 
ity for the 50.10 case with h = M/40 resolution. The 
mode-overlaps have been sorted from largest to smallest. 
The total recoil is labeled with an inverted triangle. Pos- 
itive mode-overlap contributions are labeled with circles 
and negative with diamonds. There are two important 
points to take from this figure: A) Note how quickly the 
contribution to V x and V v from each mode-overlap falls 
off; that is, there are few mode-overlaps that have signif- 
icant contribution. B) The two most dominant mode- 
overlaps (2,— 2|2, — 1) and (2,2|2, 1) contribute almost 
equally 54% (note that other modes contribute nega- 
tively) in V x and 40 % in V y . 
Another way of showing the dominance of the 
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Fig. 6. — Contribution to the recoil velocity components V x 
and V y from each (£,m\£, fh) mode-overlap for the 50.10 case with 
resolution h = M/40 extracted at = 30 M. The recoil from 
combining all mode-overlaps is labeled with an inverted triangle. 
Positive mode-overlap contributions are labeled with circles and 
negative with diamonds. 




150 
time (M) 

Fig. 7. — Recoil velocity components V x and V y versus time 
for the case SO. 10 case with resolution h = M/40 extracted at 
r det = 30 M. The solid line gives the accumulation in time of re- 
coil from all mode-overlaps combined, dotted line denotes the com- 
bined accumulations of only the two most dominant mode-overlaps, 
(2, — 2|2, —1) and (2, 2|2, 1), and the dashed line the accumulation 
in time of the (2, — 2|2, —1) mode-overlap. 

(2, — 2|2, — 1) and (2,2|2,1) mode-overlaps is presented 
in Figure [7] This figure shows the accumulated veloc- 
ity as a function of time. The solid line gives the accu- 
mulation in time of recoil from all mode-overlaps com- 
bined, the dotted line shows the combined accumulations 
of the two most dominant mode-overlaps, (2, — 2|2, — 1) 
and (2, 2 1 2, 1), and the dashed line displays the accumu- 
lation in time of the (2, — 2|2, — 1) mode overlap. 

Given that the mode-overlaps (2, — 2|2, — 1} and 
(2,2|2,1) are the principal contributors to the total 
kick velocity, we analyzed the convergence properties 
of these overlaps. Figure [5] displays the differences of 
the (2, — 2|2,— 1) mode-overlap from three resolutions, 
h/M = (1/32,1/35,1/40). The solid line is the differ- 
ence between the coarse and medium resolutions ("c- 
m"). The other lines show the difference between the 
medium and fine resolutions ("m-f"), scaled to match 
("c-m") for 3rd, 4th and 5th order convergence. It is 
clear from this figure that this mode-overlap is close to 
being 4th-order convergent. A similar situation occurs 
for the other dominant mode-overlap (2, 2|2, 1). Unfor- 
tunately, the situation is different for the other weaker 
mode-overlaps. These overlaps involve higher modes of 
\l/4 that are much more difficult to resolve given the range 
of resolutions we have. When these weaker modes are 
added to obtain the total recoil, one is no longer able 
to reach the desired 4th-order convergence. In some in- 
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Fig. 8. — Convergence analysis of the recoil contribution from 
the dominant overlap (2, — 2|2, — 1) for the SO. 10 case extracted 
a t r<i e t = 30 M. The solid line gives the difference between the 
coarse and medium resolutions ("c-m"). The other lines show the 
difference between the medium and fine resolutions ( "m-f " ) , scaled 
to match ("c-m") for 3rd, 4th and 5th order convergence. 

stances it drops to lst-order convergence. Fortunately, 
as we have seen from Figure [6j their contribution to the 
overall recoil is small. We are confident that our total 
kick velocities will not change significantly if one is able 
to achieve finer resolutions than h = M/40. 

6. SUMMARY AND DISCUSSION 

We have computed estimates of BH merger kick ve- 
locities from previously untreated physical effects arising 
from the spin of the holes. Our computational simula- 
tions provided firm predictions of kick velocities for BBH 
systems of equal mass and anti-aligned spins. Because 
we are able to accurately resolve the dominant modes 
that contribute to the kick and estimate those kicks by 
a number of methods, we are confident in our astrophys- 
ical conclusions involving the binary types we consid- 
ered. Previous studies which considered the merger of 
(non-spinning) BHs of unequal masses produced kicks 

4 Soon after the completion of our work, result s that support our 
findings of spin effects on kicks were obtain by Campanclli ct al. 
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~ 200 km s 1 with a reasonably broad maximum near 
the symmetrized mass ratio of r\ = 0.2 (mass ratio 0.38). 
From the astrophysical point of view, 200 km s~ is inter- 
esting. For instance, the escape velocity from the center 
of dwarf elliptical galaxies is 300 km s _1 , assuming the 
standard picture of dark matter halos. We found spin 
kick velocities V = 475 km s 1 a , where a is the dimen- 
sionless spin parameter, in opposite-spin configurations 
(see Figure [SJ. 

For black holes (10 — 20 M©) seen in the galaxy, 
there a re observations suppo rting spin parameters 
a £ 0.8 (McCl intock et afl 12006). and theoretical expla- 
nations of why this is so are generally applicable to 
SMBHs also. Thus we expect substantial kicks due to 
spin interactions. Our simulations predict typical kicks 
zi, 400 km s^ 1 in astrophysical BH mergers of all masses. 
These results could explain the observed absence of cen- 
tral black holes in dwarf elliptical galaxies. Our simu- 
lations show limitations, mostly due to the high cost of 
performing very high resolution runs. But, because we 
are able to accurately resolve the dominant modes that 
contribute to the kick, we believe that our astrophysical 
conclusions are secure . 
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